#ifndef __DIVTP13D3D__
#define __DIVTP13D3D__
#include <stdbool.h>
#include "mesh.h"
#include "enumerations.h"
#include "constants.h"
// ***********************************************************************
// Hdiv P1 element, DIV CONFORMING EDGE ELEMENTS OF NÉDÉLEC, 3D
// ***********************************************************************
static INT Div_T_P1_3D_3D_dof[4] = {0, 0, 1, 0};
static INT Div_T_P1_3D_3D_Num_Bas = 4;
static INT Div_T_P1_3D_3D_Value_Dim = 3;
static INT Div_T_P1_3D_3D_Inter_Dim = 1;
static INT Div_T_P1_3D_3D_Polydeg = 1;
static bool Div_T_P1_3D_3D_IsOnlyDependOnRefCoord = 0;
static INT Div_T_P1_3D_3D_Accuracy = 1;
static MAPTYPE Div_T_P1_3D_3D_Maptype = Div;

static void Div_T_P1_3D_3D_InterFun(DOUBLE RefCoord[], INT dim, DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 2 * x;
   values[1] = 2 * y;
   values[2] = 2 * z;
   values[3] = 2 * x - 2;
   values[4] = 2 * y;
   values[5] = 2 * z;
   values[6] = 2 * x;
   values[7] = 2 * y - 2;
   values[8] = 2 * z;
   values[9] = 2 * x;
   values[10] = 2 * y;
   values[11] = 2 * z - 2;
}
static void Div_T_P1_3D_3D_D000(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 2 * x;
   values[1] = 2 * y;
   values[2] = 2 * z;
   values[3] = 2 * x - 2;
   values[4] = 2 * y;
   values[5] = 2 * z;
   values[6] = 2 * x;
   values[7] = 2 * y - 2;
   values[8] = 2 * z;
   values[9] = 2 * x;
   values[10] = 2 * y;
   values[11] = 2 * z - 2;
}
static void Div_T_P1_3D_3D_D100(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 2;
   values[1] = 0;
   values[2] = 0;
   values[3] = 2;
   values[4] = 0;
   values[5] = 0;
   values[6] = 2;
   values[7] = 0;
   values[8] = 0;
   values[9] = 2;
   values[10] = 0;
   values[11] = 0;
}
static void Div_T_P1_3D_3D_D010(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 0;
   values[1] = 2;
   values[2] = 0;
   values[3] = 0;
   values[4] = 2;
   values[5] = 0;
   values[6] = 0;
   values[7] = 2;
   values[8] = 0;
   values[9] = 0;
   values[10] = 2;
   values[11] = 0;
}
static void Div_T_P1_3D_3D_D001(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 0;
   values[1] = 0;
   values[2] = 2;
   values[3] = 0;
   values[4] = 0;
   values[5] = 2;
   values[6] = 0;
   values[7] = 0;
   values[8] = 2;
   values[9] = 0;
   values[10] = 0;
   values[11] = 2;
}
static void Div_T_P1_3D_3D_D200(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 0;
   values[1] = 0;
   values[2] = 0;
   values[3] = 0;
   values[4] = 0;
   values[5] = 0;
   values[6] = 0;
   values[7] = 0;
   values[8] = 0;
   values[9] = 0;
   values[10] = 0;
   values[11] = 0;
}
static void Div_T_P1_3D_3D_D020(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 0;
   values[1] = 0;
   values[2] = 0;
   values[3] = 0;
   values[4] = 0;
   values[5] = 0;
   values[6] = 0;
   values[7] = 0;
   values[8] = 0;
   values[9] = 0;
   values[10] = 0;
   values[11] = 0;
}
static void Div_T_P1_3D_3D_D002(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 0;
   values[1] = 0;
   values[2] = 0;
   values[3] = 0;
   values[4] = 0;
   values[5] = 0;
   values[6] = 0;
   values[7] = 0;
   values[8] = 0;
   values[9] = 0;
   values[10] = 0;
   values[11] = 0;
}
static void Div_T_P1_3D_3D_D110(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 0;
   values[1] = 0;
   values[2] = 0;
   values[3] = 0;
   values[4] = 0;
   values[5] = 0;
   values[6] = 0;
   values[7] = 0;
   values[8] = 0;
   values[9] = 0;
   values[10] = 0;
   values[11] = 0;
}
static void Div_T_P1_3D_3D_D101(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 0;
   values[1] = 0;
   values[2] = 0;
   values[3] = 0;
   values[4] = 0;
   values[5] = 0;
   values[6] = 0;
   values[7] = 0;
   values[8] = 0;
   values[9] = 0;
   values[10] = 0;
   values[11] = 0;
}
static void Div_T_P1_3D_3D_D011(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x = RefCoord[0];
   y = RefCoord[1];
   z = RefCoord[2];
   values[0] = 0;
   values[1] = 0;
   values[2] = 0;
   values[3] = 0;
   values[4] = 0;
   values[5] = 0;
   values[6] = 0;
   values[7] = 0;
   values[8] = 0;
   values[9] = 0;
   values[10] = 0;
   values[11] = 0;
}
static void Div_T_P1_3D_3D_Nodal(ELEMENT *elem, FUNCTIONVEC *fun, INT dim, DOUBLE *values)
{
   // 直接给出积分格式
   INT NumPoints = 13;
   DOUBLE QuadX213[13] = {0.333333333333000,
                          0.260345966079000,
                          0.260345966079000,
                          0.479308067842000,
                          0.065130102902000,
                          0.065130102902000,
                          0.869739794196000,
                          0.312865496005000,
                          0.638444188570000,
                          0.048690315425000,
                          0.638444188570000,
                          0.312865496005000,
                          0.048690315425000};
   DOUBLE QuadY213[13] = {0.333333333333000,
                          0.260345966079000,
                          0.479308067842000,
                          0.260345966079000,
                          0.065130102902000,
                          0.869739794196000,
                          0.065130102902000,
                          0.638444188570000,
                          0.048690315425000,
                          0.312865496005000,
                          0.312865496005000,
                          0.048690315425000,
                          0.638444188570000};
   DOUBLE QuadW213[13] = {-0.149570044468000,
                          0.175615257433000,
                          0.175615257433000,
                          0.175615257433000,
                          0.053347235609000,
                          0.053347235609000,
                          0.053347235609000,
                          0.077113760890000,
                          0.077113760890000,
                          0.077113760890000,
                          0.077113760890000,
                          0.077113760890000,
                          0.077113760890000};
   INT i, j, k;
   INT num_fun = dim / 3;
   // 这里的法向量是利用$line(ii0->ii1)\times line(ii0->ii2)$生成的（为了保证同一个面在不同单元中自由度是一致的）
   INT ii0[4] = {1, 0, 0, 0};
   INT ii1[4] = {2, 2, 1, 1};
   INT ii2[4] = {3, 3, 3, 2};
   INT ii3[4] = {0, 1, 2, 3};
   INT vert04face[4] = {1, 0, 0, 0};
   INT vert14face[4] = {2, 2, 1, 1};
   INT vert24face[4] = {3, 3, 3, 2};
   DOUBLE vec_line0[3], vec_line1[3], vec_norm[3]; // 面的法向量(按照单元内部的方向)
   DOUBLE *values_temp = malloc(dim * sizeof(DOUBLE));
   DOUBLE coord[3];
   for (i = 0; i < 4; i++)
   {
      for (k = 0; k < num_fun; k++)
      {
         values[k + i * num_fun] = 0.0;
      }
      // 计算面的法向量(并非单位法向量 而是$line(ii0->ii1)\times line(ii0->ii2)$ 向量长度 = 对应面的面积*2)
      vec_line0[0] = elem->Vert_X[ii1[i]] - elem->Vert_X[ii0[i]];
      vec_line0[1] = elem->Vert_Y[ii1[i]] - elem->Vert_Y[ii0[i]];
      vec_line0[2] = elem->Vert_Z[ii1[i]] - elem->Vert_Z[ii0[i]];
      vec_line1[0] = elem->Vert_X[ii2[i]] - elem->Vert_X[ii0[i]];
      vec_line1[1] = elem->Vert_Y[ii2[i]] - elem->Vert_Y[ii0[i]];
      vec_line1[2] = elem->Vert_Z[ii2[i]] - elem->Vert_Z[ii0[i]];
      vec_norm[0] = vec_line0[1] * vec_line1[2] - vec_line0[2] * vec_line1[1];
      vec_norm[1] = vec_line0[2] * vec_line1[0] - vec_line0[0] * vec_line1[2];
      vec_norm[2] = vec_line0[0] * vec_line1[1] - vec_line0[1] * vec_line1[0];
      for (j = 0; j < NumPoints; j++)
      {
         // coord
         coord[0] = (1.0 - QuadX213[j] - QuadY213[j]) * elem->Vert_X[vert04face[i]] + QuadX213[j] * elem->Vert_X[vert14face[i]] + QuadY213[j] * elem->Vert_X[vert24face[i]];
         coord[1] = (1.0 - QuadX213[j] - QuadY213[j]) * elem->Vert_Y[vert04face[i]] + QuadX213[j] * elem->Vert_Y[vert14face[i]] + QuadY213[j] * elem->Vert_Y[vert24face[i]];
         coord[2] = (1.0 - QuadX213[j] - QuadY213[j]) * elem->Vert_Z[vert04face[i]] + QuadX213[j] * elem->Vert_Z[vert14face[i]] + QuadY213[j] * elem->Vert_Z[vert24face[i]];
         // u(t_j)
         fun(coord, dim, values_temp);
         // 面上有一个自由度
         // \int_{f} \vec{u}\cdot \vec{n}df = \sum_j \vec{u}(t_j)cdot \vec{n}\omega_j|f| = \sum_j \vec{u}(t_j)cdot (|f|\vec{n}) \omega_j
         for (k = 0; k < num_fun; k++)
         {
            values[k + i * num_fun] += (values_temp[3 * k] * vec_norm[0] + values_temp[3 * k + 1] * vec_norm[1] + values_temp[3 * k + 2] * vec_norm[2]) * QuadW213[j] / 2.0;
         }
      }
   }
   // 用面3来判断det的正负 混合积 （01\times02）\cdot 03
   if ((vec_norm[0] * (elem->Vert_X[ii3[i]] - elem->Vert_X[ii0[i]]) + vec_norm[1] * (elem->Vert_Y[ii3[i]] - elem->Vert_Y[ii0[i]]) + vec_norm[2] * (elem->Vert_Z[ii3[i]] - elem->Vert_Z[ii0[i]])) > 0)
   {
      elem->is_det_positive = 1;
   }
   else
   {
      elem->is_det_positive = 0;
   }
   free(values_temp);
}
#endif
